Methods and systems for hex-mesh optimization via edge-cone rectification

ABSTRACT

Methods and systems improve quality of a hex-mesh by: performing a first iterative procedure which starts with the input hex-mesh and which outputs an output hex-mesh having improved quality. Performing the first iterative procedure comprises: initializing a current hex-mesh to be equal to the input hex mesh; for each iteration of the first iterative procedure: performing an optimization of an energy function over a plurality of directed-edges in the current hex-mesh to determine updated vertex positions for vertices in the current hex-mesh, wherein for each directed edge, the energy function comprises a term that expresses a preference for the directed edge to be aligned with normal vectors of base triangles of an edge-cone corresponding to the directed edge; and updating the current hex-mesh with the updated vertex positions; and after one or more iterations, setting the output hex-mesh to be equal to the current hex-mesh.

REFERENCE TO RELATED APPLICATIONS

This application claims the benefit of the priority, and the benefit of35 USC 119(e), of U.S. application Ser. No. 62/196894 which is herebyincorporated herein by reference.

TECHNICAL FIELD

The invention relates to computer-based hexahedral mesh (hex-mesh)modeling of three-dimensional objects. Particular embodiments providemethods and systems for optimizing hex-meshes via edge-conerectification.

BACKGROUND

Many computer-based engineering models use finite element discretizationbased models to model and simulate the behavior of three-dimensionalobjects. A popular choice for many such engineering models involves theuse of hexahedral meshes (hex-meshes) as the finite elementdiscretization of choice. Typical, but non-limiting examples of inputobjects which may be modelled using hex-meshes include: parts ofmechanical systems; fluid conduits for fluid simulation; biologicalelements, such as human bones and/or organs. There is a general desireto provide high quality hex-meshes, since the quality of such hex-meshescan influence the accuracy and/or reliability of the correspondingengineering models. The generation of hex-meshes is typically a two-stepprocess which involves: (i) generating an initial hex-mesh whoseconnectivity is designed to fit the input object; and (ii) modifying thepositions of vertices to optimize the shapes of the mesh elements whilekeeping the connectivity fixed.

A challenge that is faced by automated (computer-based) hex-meshgeneration techniques is that the generated hex-meshes typically containpoorly shaped hexahedrons. For example, poorly shaped hexahedrons maycomprise hexahedrons that are relatively far from cuboid in shape and/orhexahedrons that are inverted or concave in shape. Even a singleinverted (concave) hexahedral element can render an entire hex-meshrepresentation unusable for simulation. Hex-meshes comprising suchpoorly shaped hexahedrons may be referred to as low quality hex-meshesand individual poorly shaped hexahedrons within such hex-meshes may bereferred to as low quality hexahedrons or low quality elements. There isa general desire for automated computer-based techniques for improvingthe quality of hex-meshes (e.g. hex-mesh representations of inputobjects). Such hex-mesh improvement techniques may improve the averagequality of the hex-mesh elements and/or the minimum quality of thehex-mesh elements and may be referred to as hex-mesh optimizationtechniques.

The foregoing examples of the related art and limitations relatedthereto are intended to be illustrative and not exclusive. Otherlimitations of the related art will become apparent to those of skill inthe art upon a reading of the description and a study of the drawings.

BRIEF DESCRIPTION OF DRAWINGS

Exemplary embodiments are illustrated in referenced figures of thedrawings. It is intended that the embodiments and figures disclosedherein are to be considered illustrative rather than restrictive.

FIG. 1 schematically illustrates an automated (computer-implemented)method for optimizing (i.e. improving the quality of) an input hex-meshto generate an improved quality output hex-mesh according to aparticular embodiment.

FIG. 2A shows an exemplary hex-segment of a hex-mesh and itscorresponding corner tetrahedrals. FIG. 2B shows an exemplary directededge in a hex-mesh and the hexahedrons surrounding the directed edge.FIG. 2C shows an exemplary directed edge in a hex-mesh, itscorresponding edge-cone, the hexahedra and corresponding base trianglesof its corresponding edge-cone and the normal vectors of its basetriangles. FIG. 2D schematically depicts a directed edge e_(ij), thebase triangles t₀, t₁, t₂, t₃ of its edge-cone and the normal vectors ofthe base triangles.

FIG. 3 is a schematic illustration of an automated(computer-implemented) method which may be used to implement the FIG. 1optimization procedure according to an example embodiment.

FIG. 4 schematically illustrates an automated (computer-implemented)method for optimizing (i.e. improving the quality of) an input hex-meshto generate an improved quality output hex-mesh according to aparticular embodiment.

FIG. 5 is a schematic illustration of an automated(computer-implemented) method which may be used to implement the FIG. 4local procedure according to an example embodiment.

FIG. 6 schematically illustrates an automated (computer-implemented)method which may be used to implement the FIG. 4 global procedureaccording to an example embodiment.

FIG. 7A shows an example of a number of topologically parallel edges ofa hexahedron. FIG. 7B shows a number of examples of topologicallyconsecutive edges of neighboring hexahedra.

FIG. 8 is a schematic depiction of a system which may be used to performany of the methods described herein according to a particularembodiment.

SUMMARY

The following embodiments and aspects thereof are described andillustrated in conjunction with systems, tools and methods which aremeant to be exemplary and illustrative, not limiting in scope. Invarious embodiments, one or more of the above-described problems havebeen reduced or eliminated, while other embodiments are directed toother improvements.

The invention provides a number of non-limiting aspects. Non-limitingaspects of the invention comprise the following:

-   1. A method, performed by a computer, for improving quality of a    hex-mesh, the method comprising:

obtaining an input hex-mesh;

performing a first iterative procedure which starts with the inputhex-mesh and which outputs an output hex-mesh having improved qualityrelative to the input hex-mesh;

wherein performing the first iterative procedure comprises:

-   -   initializing a current hex-mesh to be equal to the input hex        mesh;    -   for each iteration of the first iterative procedure:        -   performing an optimization (e g minimization) of an energy            function over a plurality of directed-edges in the current            hex-mesh to determine updated vertex positions for vertices            in the current hex-mesh, wherein for each directed edge, the            energy function comprises a term that expresses a preference            for the directed edge to be aligned with normal vectors of            base triangles of an edge-cone corresponding to the directed            edge; and        -   updating the current hex-mesh with the updated vertex    -   positions; and after one or more iterations, setting the output        hex-mesh to be equal to the current hex-mesh.

-   2. A method according to aspect 1 or any other aspect herein    wherein, for each directed edge extending from a base vertex v_(i)    to an end vertex v_(j):

the edge-cone corresponding to the directed edge comprises a pluralityof tetrahedra surrounding the directed edge, each of the plurality oftetrahedra bounded by the directed edge and by a pair of distinct edgesof a hexahedron in the current hex-mesh mesh, each of the pair ofdistinct edges emanating from the base vertex v_(i) and not includingthe end vertex v_(j);

the base triangles of the edge-cone corresponding to the directed edgecomprise, for each tetrahedron of the plurality of tetrahedra, thetriangular surface of the tetrahedron which includes the base vertexv_(i) but does not include the end vertex v_(j).

-   3. A method according to any one of aspects 1 and 2 or any other    aspect herein wherein, for each directed edge, the term of the    energy function that expresses the preference for the directed edge    to be aligned with normal vectors of base triangles of the edge-cone    corresponding to the directed edge has the form of

${\sum\limits_{t_{k} \in {T{(e_{ij})}}}\left( {\frac{e_{ij}}{e_{ij}} - n_{k}} \right)^{2}},$

as described above in connection with equation (2).

-   4. A method according to any one of aspects 1 to 3 or any other    aspect herein comprising, for each directed edge, at least one of:

determining the normal vectors of the base triangles corresponding tothe directed edge according to an equation having the form of equation(4) described above; and

obtaining target edge lengths for the output hex-mesh and determiningthe normal vectors of the base triangles corresponding to the directededge according to an equation having the form of equation (6) describedabove.

-   5. A method according to aspect 4 or any other aspect herein wherein    obtaining target edge lengths for the output hex-mesh comprises one    or more of: receiving the target edge lengths for the output    hex-mesh; using the edge-lengths of the input hex-mesh as the target    edge lengths for the output hex-mesh; and estimating the target edge    lengths for the output hex-mesh.-   6. A method according to aspect 5 or any other aspect herein wherein    obtaining target edge lengths for the output hex-mesh comprises    estimating the target edge lengths for the output hex-mesh and    estimating the target edge lengths for the output hex-mesh comprises    performing an edge-length optimization which optimizes (e.g.    minimizes) an edge-length energy function over a plurality of edges    in either the input hex-mesh or the current hex-mesh.-   7. A method according to aspect 6 or any other aspect herein wherein    the edge-length energy function comprises a first term which    expresses a preference for at least some of the plurality of edges    to have target edge lengths that are similar to the average edge    length of the plurality of edges.-   8. A method according to any one of aspects 6 to 7 or any other    aspect herein where the edge-length energy function comprises a    surface term which expresses a preference for at least some of the    plurality of edges on a boundary surface to have target edge lengths    that are similar to their original edge lengths.-   9. A method according to any one of aspects 6 to 8 or any other    aspect herein where the edge-length energy function comprises a    topologically parallel term which expresses a preference for at    least some topologically parallel edges to have target edge lengths    that are similar to one another.-   10. A method according to any one of aspects 6 to 9 or any other    aspect herein where the edge-length energy function comprises a    topologically consecutive term which expresses a preference for at    least some topologically consecutive edges to have target edge    lengths that are similar to one another.-   11. A method according to any one of aspects 6 to 10 or any other    aspect herein where the edge-length energy function has the form of    equation (22) described herein.-   12. A method according to any one of aspects 1 to 11 or any other    aspect herein wherein performing the optimization of the energy    function is subject to non-inversion constraints which require at    least one of: the hexahedra in the current hex-mesh to be    non-inverted and the tetrahedra corresponding to each directed edge    to be non-inverted.-   13. A method according to aspect 12 or any other aspect herein    wherein the non-inversion constraints have the form of equation (3)    described herein.-   14. A method according to any one of aspects 1 to 11 or any other    aspect herein wherein performing the optimization of the energy    function comprises performing an iterative optimization and wherein,    for each iteration of the iterative optimization, the energy    function comprises, for each directed edge, a weighted penalty term    which expresses a preference for the directed edge to be aligned    with a target edge direction {circumflex over (n)}_(ij).-   15. A method according to aspect 14 or any other aspect herein    wherein, for each iteration of the iterative optimization, updating    the weighted penalty term for each directed edge.-   16. A method according to any one of aspect 14 to 15 or any other    aspect herein wherein performing the first iterative procedure    comprises, for each iteration of the first iterative procedure and    for each directed edge, determining the target edge direction    {circumflex over (n)}_(ij).-   17. A method according to aspect 16 or any other aspect herein    wherein, for each directed edge, determining the target edge    direction {circumflex over (n)}_(ij) is based on characteristics of    the current hex-mesh local to the edge-cone corresponding to the    directed edge.-   18. A method according to aspect 17 or any other aspect herein    wherein, for each directed edge, determining the target edge    direction {circumflex over (n)}_(ij) is based only on    characteristics of the current hex-mesh local to the edge-cone    corresponding to the directed edge.-   19. A method according to any one of aspects 17 to 18 or any other    aspect herein wherein, for each directed edge, the characteristics    of the current hex-mesh local to the edge-cone corresponding to the    directed edge are the normal vectors of the base triangles of the    edge-cone corresponding to the directed edge and the positions of    the vertices which define the directed edge.-   20. A method according to any one of aspects 16 to 19 or any other    aspect herein wherein, for each directed edge, determining the    target edge direction {circumflex over (n)}_(ij) comprises    performing a local optimization which optimizes (e g minimizes) a    local energy function comprising a plurality of local terms, each    local term expressing a preference for the target edge direction    {circumflex over (n)}_(ij) to be aligned with a normal vector of a    corresponding one of the base triangles of the edge-cone    corresponding to the directed edge.-   21. A method according to aspect 20 or any other aspect herein    wherein, for each directed edge, performing the local optimization    comprises determining the target edge direction {circumflex over    (n)}_(ij) in accordance with equation having the form of (9)    described above.-   22. A method according to any one of aspects 20 to 21 or any other    aspect herein wherein, for each directed edge, the local    optimization is subject to non-inversion constraints.-   23. A method according to aspect 22 or any other aspect herein    wherein, for each directed edge, the non-inversion constraints have    the form of at least one of: equation (3) described herein;    equations (10) and (11) described herein; or equations (26) and (11)    described herein.-   24. A method according to any one of aspects 16 to 19 or any other    aspect herein wherein, for each directed edge, determining the    target edge direction {circumflex over (n)}_(ij) comprises:

initially setting an initial target edge direction {circumflex over(n)}_(ij) to correspond to the direction of an average of the normalvectors of the base triangles of the edge-cone corresponding to thedirected edge;

evaluating whether the initial target edge direction {circumflex over(n)}_(ij) satisfies non-inversion constraints; and

if the initial target edge direction {circumflex over (n)}_(ij)satisfies non-inversion constraints, then determining the target edgedirection {circumflex over (n)}_(ij) to be the initial target edgedirection {circumflex over (n)}_(ij) and, otherwise, performing a localoptimization which minimizes a local energy function which comprising aplurality of local terms, each local term expressing a preference forthe target edge direction {circumflex over (n)}_(ij) to be aligned witha normal vector of a corresponding one of the base triangles of theedge-cone corresponding to the directed edge.

-   25. A method according to aspect 24 or any other aspect herein    wherein, for each directed edge, performing the local optimization    comprises determining the target edge direction {circumflex over    (n)}_(ij) in accordance with equation (9) described above.-   26. A method according to any one of aspects 24 to 25 or any other    aspect herein wherein, for each directed edge, the local    optimization is subject to the non-inversion constraints.-   27. A method according to any one of aspects 24 to 26 or any other    aspect herein wherein, for each directed edge, the non-inversion    constraints have the form of at least one of: equation (3) described    herein; equations (10) and (11) described herein; or equations (26)    and (11) described herein.-   28. A method according to any one of aspects 16 to 27 or any other    aspect herein wherein, for each iteration of the iterative    optimization, the weighted penalty term is based on the target edge    direction {circumflex over (n)}_(ij).-   29. A method according to aspect 28 or any other aspect herein    wherein, for each directed edge, the weighted penalty term comprises    a penalty function E_(ij) ^(p) which is optimized as part of the    energy function and a penalty weight w_(ij) ^(p) that is updated    between iterations of the second iterative optimization.-   30. A method according to aspect 29 or any other aspect herein    wherein, for each directed edge, the penalty function E_(ij) ^(p)    expresses the preference for the directed edge to be aligned with a    target edge direction {circumflex over (n)}_(ij).-   31. A method according to any one of aspects 29 to 30 or any other    aspect herein wherein, for each directed edge, the penalty function    E_(ij) ^(p) is a quadratic function of the target edge direction    {circumflex over (n)}_(ij).-   32. A method according to any one of aspects 29 to 31 or any other    aspect herein wherein, for each directed edge, the penalty function    E_(ij) ^(p) has the form of equation (14) described above.-   33. A method according to any one of aspects 29 to 32 or any other    aspect herein wherein, for each directed edge, the penalty weight    w_(ij) ^(p) does not decrease between iterations of the iterative    optimization.-   34. A method according to any one of aspects 29 to 33 or any other    aspect herein wherein, for each directed edge, the penalty weight    w_(ij) ^(p) is relatively low when the directed edge is relatively    more aligned with a target edge direction {circumflex over (n)}_(ij)    and relatively high when the directed edge is relatively less    aligned with the edge direction {circumflex over (n)}_(ij).-   35. A method according to any one of aspects 29 to 34 or any other    aspect herein wherein the penalty weight w_(ij) ^(p) has the form of    equation (16) described above.-   36. A method according to any one of aspects 28 to 35 or any other    aspect herein wherein, for each directed edge, the weighted penalty    term comprises a geometry factor W(α_(ij)) based at least in part on    a largest angle α_(ij) between the target direction {circumflex over    (n)}_(ij) and the normal vectors of base triangles of the edge-cone    corresponding to the directed edge.-   37. A method according to aspect 36 or any other aspect herein    wherein the geometry factor W(a_(u)) increases as the largest angle    a_(u) grows.-   38. A method according to any one of aspects 36 to 37 or any other    aspect herein wherein the geometry factor W(a_(u)) has the form of    either: equation (17) described above; or equation (27) described    above.-   39. A method according to any one of aspects 14 to 38 or any other    aspect herein wherein, for each iteration of the iterative    optimization, the energy function comprises a parallel-edges    regularization term E_(regularize) _(_) _(parallel), which expresses    a preference for directions of topologically parallel edges of the    same hexahedron to be the same.-   40. A method according to aspect 39 or any other aspect herein    wherein the parallel-edges regularization term E_(regularize) _(_)    _(parallel), has the form of the left-most expression on the    right-hand side of equation (19) described herein.-   41. A method according to any one of aspects 14 to 40 or any other    aspect herein wherein, for each iteration of the iterative    optimization, the energy function comprises a consecutive-edges    regularization term E_(regularize) _(_) _(consecutive), which    expresses a preference for directions of topologically consecutive    edges of face-adjacent hexahedra to be the same.-   42. A method according to aspect 41 or any other aspect herein    wherein the consecutive-edges regularization term E_(regularize)    _(_) _(consecutive) has the form of the right-most expression on the    right-hand side of equation (19) described herein.-   43. A method according to any one of aspects 20 to 27 or any other    aspect herein wherein, for each surface directed edge having both of    its vertices on a surface boundary of the input mesh, performing the    local optimization comprises optimizing (e.g. minimizing) a    boundary-edge energy function which expresses a preference for the    surface directed edge to remain on the surface boundary.-   44. A method according to aspect 43 or any other aspect herein    wherein optimizing the boundary-edge energy function comprises    solving an equation having the form of equation (20) described    herein.-   45. A method according to any one of aspects 1 to 44 or any other    aspect herein wherein the energy function comprises a boundary    preservation term E_(boundary) which expresses a preference for    vertices which are on a surface boundary of the input hex mesh to    remain on the boundary surface.-   46. A method according to any one of aspects 14 to 44 or any other    aspect herein wherein, for each iteration of the iterative    optimization, the energy function comprises a boundary preservation    term E_(boundary) which expresses a preference for vertices which    are on a surface boundary of the input hex mesh to remain on the    boundary surface.-   47. A method according to any one of aspects 45 to 46 or any other    aspect herein wherein the boundary preservation term E_(boundary)    depends, for each surface vertex, whether the surface vertex is a    feature vertex of an object underlying the input mesh, a corner    vertex of the object underlying the input mesh or a regular surface    vertex of the object underlying the input mesh.-   48. A method according to aspect 47 or any other aspect herein    wherein the boundary preservation term E_(boundary) has a form of    equation (21) described above.-   49. A method according to any one of aspects 45 to 47 or any other    aspect herein wherein the boundary preservation term E_(boundary)    has a form of at least one of the terms on the right hand side of    equation (21) described above.-   50. A method according to any one of aspects 20 to 27 or any other    aspect herein wherein, for each directed edge, performing the    location optimization comprises performing the local optimization    with first non-inversion constraints until the current mesh is free    from inverted hexahedra and then switching the first non-inversion    constraints to first quality-improvement constraints which require    improved quality relative to previous iterations of the local    optimization.-   51. A method according to aspect 50 or any other aspect herein    wherein the first non-inversion constraints have the form of    equation (9) described above and the first quality improvement    constraints have the form of equation 26) described above.-   52. A method according to any one of aspect 36 to 38 or any other    aspect herein comprising, for each iteration of the first iterative    process, determining a worst angle w from among the largest angles    α_(ij) across the plurality of edge-cones corresponding to the    plurality of directed edges in the current mesh, and if w is less    than 90°, changing the geometry factor W(α_(ij)) from a first    geometry factor to a second geometry factor for subsequent    iterations of the first iterative process, the second geometry    factor different from the first geometry factor.-   53. A method according to aspect 52 or any other aspect herein    wherein the first geometry factor has the form of equation (17)    described above.-   54. A method according to any one of aspects 52 to 53 or any other    aspect herein wherein the second geometry factor has the form of    equation (27) described above.-   55. A system for improving quality of a hex-mesh, the system    comprising:

one or more processors configured to:

obtain an input hex-mesh;

perform a first iterative procedure which starts with the input hex-meshand which outputs an output hex-mesh having improved quality relative tothe input hex-mesh;

wherein the one or more processors are configured to perform the firstiterative procedure by:

-   -   initializing a current hex-mesh to be equal to the input hex        mesh;    -   for each iteration of the first iterative procedure:        -   performing an optimization (e g minimization) of an energy            function over a plurality of directed-edges in the current            hex-mesh to determine updated vertex positions for vertices            in the current hex-mesh, wherein for each directed edge, the            energy function comprises a term that expresses a preference            for the directed edge to be aligned with normal vectors of            base triangles of an edge-cone corresponding to the directed            edge; and        -   updating the current hex-mesh with the updated vertex            positions; and    -   after one or more iterations, setting the output hex-mesh to be        equal to the current hex-mesh.

-   56. A system according to aspect 55 which comprises any feature,    combination or features or sub-combination of features analogous to    any of aspects 1 to 54.

-   57. A non-transitory computer-readable medium comprising computer    executable code that, when executed by a computer system comprising    one computer or a plurality of computers operatively connected using    a data communications network, causes the computer system to perform    the method of aspect 1.

-   58. A non-transitory computer-readable medium according to aspect 57    or any other aspect herein that, when executed by a computer system    comprising one computer or a plurality of computers operatively    connected using a data communications network, causes the computer    system to perform any feature, combination or features or    sub-combination of features analogous to any of aspects 1 to 54.

In addition to the exemplary aspects and embodiments described above,further aspects and embodiments will become apparent by reference to thedrawings and by study of the following detailed descriptions.

Description

Throughout the following description specific details are set forth inorder to provide a more thorough understanding to persons skilled in theart. However, well known elements may not have been shown or describedin detail to avoid unnecessarily obscuring the disclosure. Accordingly,the description and drawings are to be regarded in an illustrative,rather than a restrictive, sense.

One aspect of the invention provides a computer-based automated methodfor improving the quality of (optimizing) an input hexahedral mesh(hex-mesh). The output of the method may be an output hex-mesh havingimproved average hexahedral element quality and/or improved minimumhexahedral element quality. As used herein, the relative quality of ahex-mesh element (i.e. a hexahedral element) may refer to how close theshape of the hex-mesh element is to a perfect cubic shape. Where a firsthex-mesh element has a shape that is relatively close to a cubic shape(as compared to the shape of a second hex-mesh element), the firsthex-mesh element may be referred to as having a higher quality (a betterquality, an improved quality and/or the like) relative to the secondhex-mesh element. Concave or inverted hex-mesh elements may beunderstood to have low quality relative to convex hex-mesh elements. Forbrevity, as used herein, the term mesh should be understood to refer toa hex-mesh, unless the context dictates otherwise.

FIG. 1 schematically illustrates an automated (computer-implemented)method 100 for optimizing (i.e. improving the quality of) an inputhex-mesh 102 to generate an output hex-mesh 158 having an improvedquality relative to input mesh 102, according to a particularembodiment. Method 100 (and its individual functional blocks) may beperformed by a suitably configured computer or computer system (referredto herein as a computer). Such a computer may comprise hardware,software, firmware or any combination thereof, and may be implementedusing any of a wide variety of components. By way of non-limitingexample, such a computer may comprise a programmed computer systemcomprising one or more processors, user input apparatus, displays and/orthe like. Such a computer may be implemented as an embedded system andmay share components (e.g. a processor) with other apparatus/devices.Such a computer may comprise one or more microprocessors, digital signalprocessors, graphics processors, field programmable gate arrays, signalconditioning circuitry and/or hardware (e.g. amplifier circuits and/orthe like) and/or the like. Components of the computer may be combined orsubdivided, and components of the computer may comprise sub-componentsshared with other components of the computer.

Method 100 commences in block 110 which comprises obtaining an inputhex-mesh 102. In general, input hex-mesh 102 can be obtained in block110 using any suitable technique. In some embodiments, block 110comprises receiving input hex-mesh 102 from an external system (e.g. viaa suitable hardware and/or software input interface, such as a USBinterface, a network interface and/or the like). In some embodiments,block 110 comprises obtaining input hex-mesh 102 from local memoryand/or from a different software application running on the samecomputer. In some embodiments, input hex-mesh 102 obtained in block 110comprises an automatically generated hex-mesh (e.g. a hex-mesh generatedby an automated (computer-based) hex-mesh generation technique).Non-limiting examples of suitable automated hex-mesh generationtechniques are disclosed, for example, in: SHEPHERD, J. F., AND JOHNSON,C. R. 2008. Hexahedral mesh generation constraints. Engineering withComputers 24, 3, 195 -213; SCHNEIDERS, R. 1996. A grid-based algorithmfor the generation of hexahedral element meshes. Engineering withComputers 12, 168 -177; MARECHAL, L. 2009. Advances in octree-basedall-hexahedral mesh generation: handling sharp features. In Proc.International Meshing Roundtable; TAUTGES, T. J., BLACKER, T., ANDMITCHELL, S. A. 1996. The whisker weaving algorithm: aconnectivity-based method for constructing all-hexahedral finite elementmeshes. Intl. Journal for Numerical Methods in Engineering 39, 19,3327-3349; MIYOSHI, K., AND BLACKER, T. 2000. Hexahedral mesh generationusing multi-axis cooper algorithm. In Proc. International MeshingRoundtable, 89-97; LI, Y., LIU, Y., XU, W., WANG, W., AND GUO, B. 2012.All-hex meshing using singularity-restricted field. ACM Trans. Graph.31, 6; GREGSON, J., SHEFFER, A., AND ZHANG, E. 2011. All-hex meshgeneration via volumetric polycube deformation. Computer Graphics Forum(Proc. SGP 2011); LIVESU, M., VINING, N., SHEFFER, A., GREGSON, J., ANDSCATENI, R. 2013. Polycut: monotone graph-cuts for polycube base-complexconstruction. ACM Trans. Graph. 32,6; and HUANG, J., JIANG, T., SHI, Z.,TONG, Y., BAO, H., AND DESBRUN, M. 2014. L₁ based construction ofpolycube maps from complex shapes. ACM Trans. Graph. 33, 3 (June),25:1-25:11.

In the illustrated embodiment of FIG. 1, block 110 also comprisesinitializing a current hex-mesh 112 to be the input hex-mesh 102.Subsequent operations in method 100 are then performed on current mesh112 which may be updated during the course of method 100.

Once input hex mesh 102 is obtained in block 110, method 100 proceeds tooptional block 120 which involves determining target edge lengths forthe hex elements of current hex-mesh 112. Optional block 120 is notnecessary. In some embodiments, block 120 is not present. Optional block120 is described in more detail below. Method 100 then proceeds to block150 which comprises a mesh optimization procedure based on an edge-coneformulation described in more detail below. The block 150 meshoptimization comprises improving the quality of input mesh 102 (or theinitialized version of current hex-mesh 112) to generate output mesh158, wherein the quality of output mesh 158 is improved relative to thatof input mesh 102.

In block 150, method 100 comprises performing an optimization procedurewhich involves determining new vertex positions for each directed edgein current hex-mesh 112 and outputting output hex-mesh 158 based onthese new vertex positions. The block 150 procedure may depend on a meshquality assessment based on evaluation of so-called edge cones. Theconcept of edge-cones and their relationship to hex-mesh quality is nowdescribed.

FIG. 2A shows an ideal (i.e. cubic) hexahedral element 200 and the eighttetrahedrals 202A-202H (collectively, tetrahedrals 202) defined by theoutgoing edges of the eight corner vertices 204A-204H (collectively,corners 204) of hex-element 200. In ideal hex-element 200, for eachcorner 204, each directed edge (not specifically enumerated) extendingfrom the corner 204 is parallel to the normal vector (or, for brevity,the normal) of the planar face that is bounded by the other edgesemanating from the corner 204. A hex-element is inverted if the anglebetween even one directed edge and its corresponding face normal (i.e.the normal of the planar face that is bounded by the other edgesemanating from the same corner) is over 90°. It can be observed fromFIG. 2A that the tetrahedrals 202 are formed by pairs of directed edgesand corresponding triangular faces, with each such pair comprising: onedirected edge that emanates from a comer; and one correspondingtriangular face (i.e. the planar triangular face which is bounded by theother edges emanating from the same corner).

FIGS. 2B and 2C schematically depict the formation of an edge-cone 212which may be used to evaluate mesh quality in accordance with particularembodiments. Each edge-cone 212 corresponds to (or has a one-to-onerelationship) with a directed edge 214. As shown best in FIG. 2C,directed edge 214 may extend from vertex v_(i) to vertex v_(j) and maybe referred to as directed edge (v_(i),v_(j)). As shown in FIG. 2B,directed edge 214 is a part of (i.e. an edge of) a number of hexahedrons210A-210E (collectively, hexahedrons 210). As shown in FIG. 2C, for eachsuch hexahedron 210, directed edge 214 may be part of (i.e. an edge of)a corresponding tetrahedron 216A-216E (collectively, tetrahedrons 216)formed by directed edge 214 and a corresponding face (i.e. the planarface which is bounded by the other edges of the same hexahedron 210which emanate from the same vertex v_(i). It should be noted that onlytetrahedrons 216A, 216B, 216E are visible in FIG. 2C. The union oftetrahedrons 216 corresponding to directed edge 214, (v_(i), v_(j)) maybe referred to as the edge-cone 212 corresponding to directed edge 214,(v_(i),v_(j)).

In some embodiments, the quality of a hex-mesh that includes a pluralityof directed edges may be expressed in terms of the differences betweenthe directions of: the directed edges which form the axes ofcorresponding edge-cones and the normal vectors of the triangular facesat their bases. For the particular edge-cone 212 shown in FIG. 2C, thequality of the edge cone may be expressed in terms of the differencesbetween the direction of directed edge directed edge 214, (v_(i), v_(j))and the normal vectors of the triangles at the bases of tetrahedrons216. By way of example, for tetrahedron 216A shown in FIG. 2C, theangular difference which goes into the quality evaluation may be theangular difference between the direction of directed edge 214, (v_(i),v_(j)) and the normal vector corresponding to the triangle formed byvertices v_(i), u_(k) and u_(k+1). FIG. 12D illustrates an exemplarydirected edge 218 and the normal vectors 220A-220D (collectively, normalvectors 220) for the triangles 222A-222D (collectively, triangles 222)at the bases of corresponding tetrahedrons (not shown in FIG. 2D). Ifthese angles are less than 90° for all the triangular faces, for all ofthe edge cones in a hex-mesh, then the hex-mesh is inversion free.

Referring back to FIG. 1, the block 150 optimization procedure maycomprise determining new vertices for each directed edge in current mesh112. Instead of directly assessing the quality of each hex-element incurrent mesh 112 and instead of assessing the quality of eachhex-element of current mesh 112 using the eight tetrahedra associatedwith the eight corners of the hex-element, the block 150 optimizationprocedure may involve evaluating the quality of current mesh 112 bytraversing the edge-cones in current mesh 112 and evaluating the qualityof each edge-cone.

A particular embodiment of block 150 is now described. For the purposesof this explanation, we define

to be a general hex-mesh with vertices

and edges

) respectively. For each directed edge (v_(i), v_(j)), we define e_(ij)to correspond to the vector v_(i)-v_(j). While the variables expressedin the description that follows may be expressed in terms of eij, insome embodiments, the variables actually optimized may comprise theunderlying vertex positions. If we consider a directed edge e_(ij) fromamong the edges

and the collection of all hexahedra H that contain e_(ij), then we maydefine a set Q={q_(k)} to be the set (e.g. the umbrella) ofquadrilaterals belonging to H and containing the vertex v_(i) but notthe vertex v_(j) we may further define T={t_(k)} to be the set (e.g. thefan) of triangles that are connected to the vertex v_(i), with theunderstanding that T has a winding order consistent with the windingorder of the quads of Q and that t_(|T|)=t₀ (i.e. the indices of t_(k)reach a finite number and then return to zero). FIG. 2D schematicallyillustrates examples of such quadrilaterals Q={q₀. . . q₃} and trianglesT={t₀ . . . t₃} for the particular example directed edge 218 shown inFIG. 2D.

Typically, an input hexahedral mesh may be assumed to have (or may beassigned) a consistent ordering of vertices, such that the vertices ofany quadrilateral forming an individual hexahedral element, when viewedfrom the center of the hexahedral element, has a counterclockwiseorientation. This in turn infers a counterclockwise orientation on thetetrahedral decomposition of the hexahedral elements of the mesh, andtherefore on the vertices of any triangle on the tetrahedron. Denotingthe vertices of an arbitrary triangle t_(k) as u_(k); v_(i); u_(k+1) (asshown for a particular example triangle in FIG. 2C), it follows that forthe tetrahedron having vertices v_(i); v_(j); u_(k); u_(k+1) to have apositive volume, the vertex v_(j) must lie on the positive side of thesupporting plane of the triangle t_(k) (i.e. the triangle havingvertices u_(k); v_(i); u_(k+1)). This positive tetrahedral volume andcorresponding condition on the vertices may be imposed for each trianglein the set T. The shape of each tetrahedron (e.g. its scaled Jacobian)is optimized if directed edge e_(ij) is aligned with the normal vectorsn_(k) corresponding to the triangles t_(k) in the set T and, when theshape of a tetrahedron is optimized, the shape of the correspondingportion of its hexahedron is also optimized.

The block 150 optimization procedure may comprise determining thedirections of each normalized directed edge ê_(ij)(e.g.

$\left. \frac{e_{ij}}{e_{ij}} \right)$

and/or corresponding positions of vertices v_(i) for current mesh 112,which minimize an energy function E_(global) comprising an energy termE_(cone), where E_(cone) depends on the difference in orientationbetween the direction of edge e_(ij) and the normal vectors n_(k) of itscorresponding triangles over each directed edge in the current mesh 112.In one particular example embodiment, the energy function E_(global)includes only the energy term E_(cone)—i.e.:

E_(global)=E_(cone)   (1a)

In another example embodiment, the energy function E_(global) is givenby:

E _(global) =E _(cone) +E _(boundary)   (1b)

where E_(boundary) represents a term that aims to preserve the boundaryof input hex-mesh 102 and which is explained in more detail below.

In one particular example embodiment, the energy term E_(cone) may begiven by:

$\begin{matrix}{E_{cone} = {\sum\limits_{e_{ij} \in E}{\sum\limits_{t_{k} \in {T{(e_{ij})}}}\left( {\frac{e_{ij}}{e_{ij}} - n_{k}} \right)^{2}}}} & (2)\end{matrix}$

where the left summation operator is over all of the directed edgese_(ij) in current hex-mesh 112 and the right summation operator is overall of the base triangles t_(k) in an edge-cone corresponding to aparticular directed edge e_(ij).

The block 150 optimization and the corresponding minimization of theenergy function E_(global) (of equation (1a) or (1b)) may be subject tonon-inversion constraints which prevent the solutions of the block 150optimization (i.e. the vertex positions v_(i)) from corresponding toinverted hex-elements and/or inverted tetrahedra. In one particularembodiment, these block 150 non-inversion constraints may have the form,for each directed edge e_(ij):

e _(ij) ·n _(k)>0 ∀e _(ij) ∈ε, t _(k) ∈T(e _(ij))   (3)

where the normal n_(k) may be explicitly defined as a function of vertexpositions:

$\begin{matrix}{{{n_{k} - \frac{\left( {u_{k} - v_{i}} \right) \times \left( {u_{k + 1} - v_{i}} \right)}{{\left( {u_{k} - v_{i}} \right) \times \left( {u_{k + 1} - v_{i}} \right)}}} = {0\; {\forall{e_{ij} \in ɛ}}}},{t_{k} \in {T\left( e_{ij} \right)}}} & (4)\end{matrix}$

It may be observed that the energy term E_(cone) of equation (2)expressly accounts for the opposite directions of each directed edgee_(ij) separately, as the tightness of the constraints on thesedirections may significantly differ.

Given the typical size of real-world meshes, optimizing the formulationof equations (1)-(4) directly using standard non-linear optimizationmethods may be impractical. In particular, the normalization by edgelength in this formulation may present optimization challenges.Consequently, in some embodiments, the block 150 optimization involvesthe use of pre-estimated edge lengths L_(ij), where L_(ij) representsthe pre-estimated length of the edge between the vertices v_(i) andv_(j). Pre-estimated edge lengths L_(ij) may be determined in block 120,which is discussed in more detail below. Substituting such pre-estimatededge lengths L_(ij) into equation (2) reduces equation (2) to thequadratic function:

$\begin{matrix}{E_{cone} = {E_{Qcone} = {\sum\limits_{e_{ij} \in E}{\sum\limits_{t_{k} \in {T{(e_{ij})}}}\left( {\frac{e_{ij}}{L_{ij}} - n_{k}} \right)^{2}}}}} & (5)\end{matrix}$

and reduces equation (4) to:

$\begin{matrix}{{{n_{k} - \frac{\left( {u_{k} - v_{i}} \right) \times \left( {u_{k + 1} - v_{i}} \right)}{{L_{k,i} \times L_{{k + 1},i}}}} = {0\; {\forall{e_{ij} \in ɛ}}}},{t_{k} \in {T\left( e_{ij} \right)}}} & (6)\end{matrix}$

FIG. 3 is a schematic illustration of a method 150A which may be used toimplement the block 150 optimization procedure according to an exampleembodiment. Method 150A comprises

-   -   (i) determining triangle normals n_(k) (in block 152) for each        triangle in the edge-cone corresponding to each directed edge        e_(ij). The block 152 triangle normal determination may be        performed using current hex-mesh 112 and may comprise using        equations having a form of either equation (4) or equation (6);    -   (ii) determining new vertex positions v₁ (in block 153) by        minimizing an energy function E _(global) having a form of        either equation (1a) or equation (1b) subject to the constraints        of equation (3). The block 153 optimization may be performed on        the basis of current hex-mesh 112. Where the approximate normals        n_(k) of equation (6) are used for step (i) (block 152), then        step (ii) (block 153) may comprise a quadratic minimization        problem with linear inequalities, solvable with standard        quadratic programming methods. The step (ii) (block 153)        quadratic minimization problem is convex and has a unique        minimum;    -   (iii) updating current hex-mesh 112 (in block 154) to have the        new vertex positions v_(i) determined in step (ii) (block 153);        and    -   (iv) repeating steps (i), (ii) and (iii) (blocks 152, 153 and        154) until satisfying a convergence criteria in block 155.        Assessing convergence in block 155 may involve evaluating one or        more convergence criteria. By way of non-limiting example, such        block 155 convergence criteria may comprise:        iteration-to-iteration changes in vertex positions that are        below a convergence threshold; iteration-to-iteration changes in        the energy function value that are below a convergence        threshold; temporal thresholds; thresholds based on a number of        iterations; until the quality of the hex mesh is above a minimum        threshold for all hexahedral elements in the hex mesh; until the        average quality of the hexahedral elements of the hex mesh is        above a certain threshold; some combination of the above; and/or        the like.    -   (v) outputting the updated current hex-mesh 112 as output hex        mesh 158 once the block 155 convergence criteria are satisfied.

With the FIG. 3 embodiment (method 150A) of the block 150 optimizationprocedure, the resultant output hex-mesh 158 at the conclusion of method150A has no inverted hex-elements. In some embodiments, when convergenceoccurs (block 155 YES branch) and the method 150A procedure terminatesby outputting the vertex positions v_(i) of current hex-mesh 112 asoutput hex-mesh 158 (block 156), the vertex positions v_(i), and hencethe corresponding triangle normal n_(k), no longer change and,accordingly, the non-inversion constraints (equation (3) hold for thenormals n_(k) of output hex-mesh 158. In some circumstances, method 150Amay terminate prematurely without finding a set of vertex positionsv_(i) that satisfies all of the non-equation (3) inversion constraintsat the same time. Such a situation can arise in a number of scenariosincluding, without limitation, the case of degenerate base triangles,whose normals computed using equation (6) have zero length. Suchscenarios can and do occur in practice if input mesh 102 is badlytangled (e.g. has a relatively large number of inverted hex-elements).However, since triangle normals can and do change as mesh vertices move,the fact that a current set of normals does not allow for a solutiondoes not necessarily imply that a solution does not exist for a givenmesh topology.

FIG. 4 schematically illustrates an automated (computer-implemented)method 300 for optimizing (i.e. improving the quality of) an inputhex-mesh 102 to generate an improved quality output hex-mesh 158according to a particular embodiment. Method 300 (and its individualfunctional blocks) may be performed by a suitably configured computerwhich may have any of the above-described characteristics of thecomputer that performs method 100. Method 300 commences in block 310which comprises obtaining an input hex-mesh 102. In general, block 310of method 300 may be substantially similar to block 110 of method 100discussed above. Like block 110 described above, block 310 may compriseinitializing a current hex-mesh 112 to be the input hex-mesh 102.Subsequent operations in method 300 are then performed on current mesh112 which may be updated during the course of method 300.

Once input hex mesh 102 is obtained and current hex-mesh 112 isinitialized in block 110, method 300 proceeds to optional block 320which involves determining target edge lengths for the hex elements ofcurrent hex-mesh 112. Optional block 320 may be substantially similar toblock 120 described elsewhere herein. In some embodiments, block 320 isnot present. Method 300 then proceeds to block 330 which comprises acombined local/global mesh optimization loop (or for brevity, a meshoptimization loop 330). Mesh optimization loop 330 comprises improvingthe quality of input hex-mesh 102 to generate output hex-mesh 158,wherein the quality of output hex-mesh 158 is improved relative to thatof input hex-mesh 102. In the illustrated FIG. 4 embodiment, meshoptimization loop 330 starts in block 360.

In block 360 of the illustrated FIG. 4 embodiment, method 300 comprisesperforming a local procedure (e.g. local to each directed edge e_(ij)and its corresponding edge-cone). Block 360 may involve determiningtarget directions {circumflex over (n)}_(ij) for each directed edgee_(ij) based only on parameters associated with the edge-conecorresponding to the directed edge e_(ij) (e.g. based on the basetriangle normals n_(k) of the edge-cone corresponding to the particulardirected edge e_(ij)). In some embodiments, block 360 may compriseperforming a local constrained optimization for each directed edgee_(ij), where the constraints comprise local non-inversion constraintsbased on the edge-cone corresponding to the particular directed edgee_(ij) (e.g. non-inversion constraints based on the base trianglenormals n_(k) of the edge cone corresponding to the particular directededge e_(ij)).

In one particular embodiment, block 360 comprises, for each directededge e_(ij), determining a target direction {circumflex over (n)}_(ij)based on a local optimization having the form:

ñ _(ij)=arg min Σ_(k=0) ^(T)(ñ _(ij) −n _(k))²   (9)

subject to the local non-inversion constraints having the form:

{circumflex over (n)} _(ij) ·n _(k)<ε<0 ∀t _(k) ∈T   (10)

∥ñ _(ij)∥²≦1   (11)

It may be observed that equations (9)-(11) form a constrainedoptimization problem whose solution {circumflex over (n)}_(ij) minimizesequation (9) subject to the constraints of equations (10) and (11). Itmay also be observed that equation (9) and (10) are only evaluated overthe local edge-cone that corresponds to the directed edge e_(ij) and thebase triangles of the local edge-cone. The convex constraint of equation(11) replaces the explicit, but non-convex, requirement for the targetdirection {circumflex over (n)}_(ij) to have unit length. In someembodiments, the equation (11) constraint could be replaced by anequality (e.g. ∥{circumflex over (n)}_(ij)∥²=1). In block 360, the basetriangle normals n_(k) may be determined based on current hex mesh 112.The current hex mesh 112 used in block 360 may be based on: theinitialized value of current hex mesh 112 (applied in block 310); or theupdate to current hex mesh 112 after the last iteration of meshoptimization loop 330 (e.g. after the block 370 global procedure, asdescribed in more detail below). The base triangle normals n_(k) may bedetermined on the basis of an equation having the form of equation (4)or equation (6) and, in some embodiments, may be normalized to have unitlength. The resulting target directions ñ_(ij) may be similarlynormalized. It should be noted that if all that was required was for thedot product {circumflex over (n)}_(ij)·n_(k) to be positive, theequation (11) constraint could be eliminated. However, in someembodiments, a greater than zero threshold ε is used in the equation(10) constraint (for reasons discussed in more detail below).Accordingly, the equation (11) constraint is explicitly included in thisdescription to prevent the resulting target directions {circumflex over(n)}_(ij) from satisfying the equation (10) constraint by scaling.

The block 360 local optimization (as characterized by the formulation ofequations (9), (10) and (11) or otherwise) may comprise, for eachdirected edge e_(ij), a quadratic convex optimization problem withlinear and quadratic inequality constraints on the three coordinates oftarget direction {circumflex over (n)}_(ij). The block 360 optimizationcan be solved efficiently and robustly using standard quadraticprogramming tools such as, by way of non-limiting example, thosedisclosed by GUROBI OPTIMIZATION, 2013 at http://www.gurobi.com/. Insome circumstances (e.g. with some directed edges of particular inputmeshes 102), equation (9) is minimized by a target direction {circumflexover (n)}_(ij) that is the conventional average of the base trianglenormal n_(k) and satisfies the non-inversion constraints (e.g.constraints of the form of equations (10) and (11)) as-is. Consequently,in some embodiments, to speed-up computation, block 360 may firstattempt to determine a base solution target direction {circumflex over(n)}_(ij) (e.g. a target direction {circumflex over (n)}_(ij) that isthe conventional average of the base triangle normal n_(k)) and evaluatewhether the base solution target direction {circumflex over (n)}_(ij)satisfies the non-inversion constraints. If it is determined that thebase solution target direction {circumflex over (n)}_(ij) does notsatisfy the non-inversion constraints, block 360 may revert to using aquadratic programming solver.

Where the block 360 local optimization is a quadratic convexoptimization problem with linear and quadratic inequality constraints, atarget direction ñ_(ij) that satisfies all the non-inversion constraintsmay not always exist. However, an advantage of a localized solution isthat method 300 can still proceed with the block 370 global procedure(or the remaining portion of an iteration of mesh optimization loop 330)in a meaningful way, even without a block 360 solution. Where, for aparticular directed edge e_(u), the block 360 local optimization doesnot yield an inversion-free solution, block 360 may, in someembodiments, return the target direction {circumflex over (n)}_(ij) tobe either the current edge direction (e.g. the direction of e_(ij) incurrent hex-mesh 112) or the average of the base triangle normals n_(k)(e.g. the unconstrained minimizer of equation (9)). In some embodiments,the decision between these two output target directions {circumflex over(n)}_(ij) may be based on the largest angle α_(ij) between the directionof the proposed output target direction {circumflex over (n)}_(ij) andthe base triangle normals n_(k). In some embodiments, where thiscircumstance arises for a particular directed edge e_(ij) in block 360,block 360 comprises selecting the output target directions {circumflexover (n)}_(ij) to be the direction for which this largest angle a_(u) isthe smallest. In some embodiments, block 360 may make this same electionof output target directions {circumflex over (n)}_(ij) if one of thebase triangle normals n_(k) has zero length. Since, in some embodiments,the block 360 local optimization seeks for all directed edges e_(ij) tohave a known target length (e.g. L_(ij)) and indirectly optimizes thebase triangle angles α_(ij), such degenerate configurations aretypically resolved during the first iteration of mesh-optimization loop330.

In some embodiments, the above discussed implementation of block 360 isused for interior edges e_(ij) (i.e. edges e_(ij) other than surfaceedge). In some embodiments, block 360 may use a modified formulation forsurface edges. This modified formulation is described in more detailbelow.

FIG. 5 schematically illustrates a method 400 which may be used toimplement the block 360 local procedure according to an exampleembodiment. Method 400 shown in FIG. 5 is implemented once for eachdirected edge e_(ij) in current hex-mesh 112 and may be based on thecharacteristics of the edge-cone corresponding to the directed edgee_(ij). Method 400 commences with block 401 which involves determiningthe triangle normals n_(k) for each base triangle of the currentdirected edge based on current hex mesh 112. The determination of thebase triangle normals n_(k) for the current directed edge in block 401may be similar to block 152 of method 150A described above and may bebased on equation (4) or equation (6) suitably modified for applicationto a single directed edge. Method 400 then proceeds to the block 402inquiry as to whether any of the base triangle normals n_(k)corresponding to the current directed edge have zero length. If thereare base triangle normals with zero length, then the block 402 inquiryis positive, and method 400 proceeds via the block 402 YES branch toblock 416 discussed further below. Otherwise, method 400 proceeds viathe block 402 NO branch to block 404, where the target direction{circumflex over (n)}_(ij) is set to be the average of the base trianglenormals (e.g. the unconstrained solution to equation (9)). Method 400then proceeds to block 406 which involves an inquiry into whether theblock 404 target direction {circumflex over (n)}_(ij) satisfies theapplicable non-inversion constraints (e.g. constraints having the formof equations (10) and (11)). If the block 406 inquiry is positive, thenthe block 404 target direction {circumflex over (n)}_(ij) is output inblock 408. If, on the other hand, the block 406 inquiry is negative,then method 400 proceeds to block 410 which involves implementing aquadratic programming solver to attempt to solve the constrainedoptimization problem (e.g. having the form of equations (9), (10) and(11)) to obtain a target direction {circumflex over (n)}_(ij). Method400 then proceeds to block 412 which involves an inquiry into whetherthe block 410 quadratic programming solver obtained a target direction{circumflex over (n)}_(ij) that satisfies the non-inversion constraints.If the block 412 inquiry is positive, then the block 420 targetdirection {circumflex over (n)}_(ij) is output in block 414. If, on theother hand, the block 412 inquiry is negative, then method 400 proceedsto block 416. Block 416 can be reached via the block 402 YES branch orthe block 412 NO branch. Block 416 may involve selecting the outputtarget direction {circumflex over (n)}_(ij) to be one of: the currentedge direction (e.g. the direction of e_(u) in current hex-mesh 112) orthe average of the base triangle normals n_(k) (e.g. the unconstrainedminimizer of equation (9) or the block 404 target direction {circumflexover (n)}_(ij)). As discussed above, in some embodiments, the particularone of the target directions {circumflex over (n)}_(ij) chosen in block416 may be based on the largest angle α_(ij) between the direction ofthe proposed output target direction {circumflex over (n)}_(ij) and thebase triangle normals n_(k). In some embodiments, block 416 comprisesselecting the output target directions {circumflex over (n)}_(ij) to bethe direction for which this angle α_(ij) is the smallest.

Referring back to FIG. 4, the block 360 local procedure outputs targetedge directions {circumflex over (n)}_(ij), 362 for each directed edgee_(ij) in current hex-mesh 112. Method 360 then proceeds to the globalprocedure in block 370. In some embodiments, the block 370 globalprocedure may comprise an iterative optimization procedure whichdetermines new vertex positions v_(i) (e.g. globally for current mesh112) that attempt to align the directed edges e_(ij) of current mesh 112with the block 360 target directions {circumflex over (n)}_(ij). Ratherthan using hard non-inversion constraints (which may result in apremature termination of the optimization procedure similar to thatdiscussed above in connection with method 100), the block 370 globalprocedure may use penalty functions between successive iterations whichmay keep the block 370 vertex positions (in each iteration) relativelyclose to the block 360 target directions {circumflex over (n)}_(ij). Theuse of penalty functions in block 370 may permit the block 370optimization to proceed to the next iteration in cases where thenon-inversion constraints cannot be immediately satisfied.

Block 370 may use an energy function E_(global) based on one ofequations (1a) or (1b) describe above and may be augmented by a penaltyterm as discussed in more detail below. In some embodiments, the block370 energy function may have the form:

E _(global) =Ê+Σ _(i) W _(i) ^(p) E _(i) ^(p)   (13)

Where: Ê is the desired energy function, but for the penalty term (e.g.the right hand side of one of equations (1a) or (1b), with the possibleaddition of other terms as described in more detail below); W_(i) ^(p)are individual penalty weights; and E_(i) ^(p) are weighted penaltyterms. Block 370 may apply per-edge-cone (or per directed edge) penaltyterms aimed to keep the edge directions determined in each iteration ofthe block 370 global procedure sufficiently close to the bock 360 targetdirections {circumflex over (n)}_(ij). In some embodiments, the block370 penalty formulation is based on the observation that the norm of thedifference between unit-length vectors is above unity when the anglebetween the unit length vectors is obtuse, and less than one, when thetwo unit-length vectors are relatively close (under 60°). Thus forlarger angles, a high order monomial of this norm may serve as a penaltywhose cost increases dramatically as the angle increases, stronglydiscouraging any solution where the new edge directions deviate farenough from the target to the point where a tetrahedron in an edge-conebecomes inverted. At the same time, the maximal distance betweenunit-length vectors is 2, providing a natural upper bound on the size ofthe penalty. With this formulation, the block 370 penalty scheme willterminate, as the penalties may not increase indefinitely.

In some embodiments, to keep the optimized energy function E_(global)quadratic, block 370 may comprise separating a penalty monomial into aquadratic term, which we treat as a penalty function

$\begin{matrix}{E_{ij}^{P} = {{\frac{e_{ij}}{{\hat{L}}_{ij}} - {\hat{n}}_{ij}}}^{2}} & (14)\end{matrix}$

and a higher order higher-order weight term w_(ij) ^(p) which is updatedbetween penalty iterations within block 370 and which, in someembodiments, is never decreased. Determination of the normalizationfactors {circumflex over (L)}_(ij) is described in more detail below. Insome embodiments, the weights w_(ij) ^(p) are initially set to 1. Aftera solution to a particular iteration within block 370 is obtained, block370 may comprise updating the penalty weights w_(ij) ^(p) by setting

$\begin{matrix}{w_{ij}^{P} = {\max \left( {w_{ij}^{P},{{\frac{e_{ij}}{{\hat{L}}_{ij}} - {\hat{n}}_{ij}}}^{6}} \right)}} & (16)\end{matrix}$

The choice of the exponent being six in equation (16) is empirical andnot mandatory. The inventors have found that lower values for thisexponent did not sufficiently penalize large angular deviations from theblock 360 desired directions {circumflex over (n)}_(ij), leading todecreases in solution quality, while larger exponents reduced numericstability (e.g. because increasing the exponent increases the ratiobetween the largest and smallest elements of the gradient matrix of theblock 370 energy function, which in turn increases its conditionnumber).

When selecting the normalization factors {circumflex over (L)}_(ij),focusing on angle only involves considering a pure directionaldifference—e.g. normalizing each directed edge of current hex-mesh 112by its length at the beginning of the block 370 global procedure. Insome circumstances, this choice may be sub-optimal (e.g. if theedge-length of a directed edge in current hex-mesh 112 is significantlydifferent from its estimated target edge length (e.g. its block 320target edge length), as the edge-length of current hex-mesh 112 may be apoor predictor of the edge-length after the block 370 optimization whichis expected to get closer to its target edge length. In someembodiments, block 370 comprises anticipating this change in edge-lengthby, for each edge, setting

${{\hat{L}}_{ij} = \frac{{e_{ij}} + L_{ij}}{2}},$

where ∥e_(ij)∥ is the edge-length of the edge in current hex-mesh 112and L_(ij) is the block 320 estimated target edge length discussedelsewhere in this description.

In some embodiments, block 370 also involves recognizing that the degreeto which an edge e_(ij) and its block 360 target direction {circumflexover (n)}_(ij) should be aligned may depend on the geometry of theedge-cones determined in the block 360 local procedure. Specifically,block 370 may comprise considering the largest angle α_(ij) between theblock 360 target direction {circumflex over (n)}_(ij) and thecorresponding edge-cone's base triangle normal n_(k) used in block 360.To minimize (or reduce) the likelihood of inversion in block 370, whenα_(ij) is below 90°, the edge e_(ij) determined in each iteration ofblock 370 and the block 360 target direction {circumflex over (n)}_(ij)should be more aligned as the angle α_(ij) grows. In some embodiments,therefore, block 370 comprises multiplying the weights w_(ij) ^(p) by afunction that penalizes deviation between the edge e_(ij) determined ineach iteration of block 370 and the block 360 target direction{circumflex over (n)}_(ij) more significantly as the angle α_(ij) grows.In some embodiment, block 370 comprises multiplying the weights w_(ij)^(p) by a function of the angle α_(ij). In one particular embodiment,block 370 comprising multiplying the weights w_(ij) ^(p) by a functionhaving the form

$\begin{matrix}{{W\left( \alpha_{ij} \right)} = {1 + {10^{\frac{- {cos\alpha}_{ij}^{2}}{2\sigma^{2}}}}}} & (17)\end{matrix}$

In some embodiments, σ is set at σ=0.3 using the three sigma rule forthe weight to drop to 1 when the angle α_(ij) is zero, allowing forlarger deviation.

In some embodiments, the block 370 global procedure may compriseiterative optimizing a quadratic energy function having the form:

E _(global) =Ê+Σ _(ij) W(α_(ij))w _(ij) ^(p) E _(ij) ^(p)   (18)

where: Ê is the desired energy function, but for the penalty term (e.g.the right hand side of one of equations (1a) or (1b), with the possibleaddition of other terms as described in more detail below);W(α_(ij))w_(ij) ^(p) are individual penalty weights; and E_(i) ^(p) areweighted penalty terms. Block 370 may comprise iteratively repeating theoptimization (e g minimizing equation (18)) with the current penaltyweights and then updating the penalty weights W (α_(ij))w_(ij) ^(p)until the vertex positions v_(i) and/or the combined energy function(e.g. equation (18)) no longer change or change less than a suitablethreshold amount. Note that, as expected from penalty terms, the weightsW(α_(ij))w_(ij) ^(p) do not decrease between iterations. It will also beappreciated that, in some embodiments, the block 370 penalties areapplied to a particular edge e_(ij) only if the block 360 targetdirections {circumflex over (n)}_(ij) satisfy the non-inversioninequalities in block 360 (corresponding to α_(ij)=90°). If thenon-inversion constraints are not satisfied for a particular edgee_(ij), then that particular edge may be omitted from penalties (e.g.omitted from the summation term of equation (18)). This can beappreciated on the basis that block 370 should only penalize changesfrom target directions corresponding to non-inverted conditions, andchanges from inverted conditions should be welcomed.

FIG. 6 schematically illustrates a method 440 which may be used toimplement the block 370 global procedure according to an exampleembodiment. Method 440 shown in FIG. 6 is implemented globally, whichmay be, but need not necessarily be, performed over the entire currentmesh 112 and/or which may be performed over a suitable portion ofcurrent mesh 112 comprising a plurality of edge-cones. Method 440commences in block 442 which comprises determining the base trianglenormal vectors n_(k) for each directed edge e_(ij) in current mesh 112.In some embodiments, these base triangle normals may already be known(e.g. from procedures performed in block 360). Method 440 then proceedsto block 444 which comprises performing an iteration of the optimizationusing the current penalty weights (e.g. W(α_(ij))w_(ij) ^(p)). Asdiscussed above, in some embodiments, the initial penalties may be setto unity on the first iteration of method 440 and may be updated in eachsuccessive iteration of method 440. The block 444 optimization maycomprise minimizing an energy function E_(global) (which may have theform of equation (18)) to determine new vertex positions v_(i). Method440 may then proceed to block 446 which comprises determining newpenalty weights (e.g. W(a_(ij))w_(ij) ^(p)). In some embodiments, theblock 446 terms w_(ij) ^(p) may be determined using the block 444 newvertex positions and an equation having the form of equation (16). Insome embodiments, the block 446 terms W (α_(ij)) may be determined usingthe block 444 new vertex positions and an equation having the form ofequation (17). Method 440 then proceeds to block 448 which involvesupdating the current hex-mesh 112 based on block 444 new vertexpositions. This block 448 updated current mesh 112 (together with theblock 446 updated weights) may be used for successive iterations ofmethod 440. Method 440 then proceeds to block 450 which involves theevaluation of one or more loop-exit (convergence) criteria. In someembodiments, such block 450 convergence criteria may comprise, withoutlimitation, evaluating whether the changes in the block 444 vertexpositions v_(i) (or equivalently changes in current mesh 112) betweensuccessive iterations are less than a suitable threshold, evaluatingwhether the changes in the combined energy function (e.g. equation (18))between successive iterations are less than a suitable threshold,evaluating a threshold based on the number of iterations, evaluating athreshold based on some temporal criteria and/or the like. If the block450 convergence criteria is not met, then method 440 loops back to block442 with the block 448 updated current hex mesh 112 and the block 446updated weights. If the block 450 convergence criteria is met, thenmethod 440 proceeds to block 452, where current hex-mesh 112 is outputas a result of method 440.

Returning back to FIG. 4, it can be seen that the current hex-mesh 112is the output of block 370 prior to an evaluation of loop exit criteriain block 380. Block 380 involves evaluation of an overall loop exitcriteria. The block 380 loop exit criteria may comprise, by way ofnon-limiting example, evaluating whether the changes in the vertexpositions v_(i) (or equivalently changes in current mesh 112) betweensuccessive iterations of mesh-optimization loop 330 are less than asuitable threshold, evaluating hard constraints (e.g. of the form ofequation (8)), evaluating a threshold based on the number of iterationsof mesh-optimization loop 330, evaluating a threshold based on sometemporal criteria and/or the like. If the block 380 loop-exit criteriaare not met, then method 300 loops again through the local/globaloptimizations of blocks 360, 370. If, however, the block 380 loop-exitcriteria are met, then method outputs current hex-mesh 112 as outputhex-mesh 158. Output hex-mesh 158 may have improved quality relative toinput hex-mesh 102.

A number of optional procedures are now described which may be used invarious embodiments and/or in various blocks of methods 100, 300.

In some embodiments, the energy function E_(global) described above maybe augmented with a regularization term E_(regularize) which may beaimed at generating more unified edge directions for adjacent mesh edgeswhich leverages the structure of the hexahedral mesh. The regularizationterm E_(regularize) may be added, for example, to the right hand side ofequation (1a), (1b) and/or (18). The regularization term E_(regularize)may form part of the term Ê in equation (18). It may be observed thatfor high quality hex-meshes, the directions of topologically paralleledges within the same hex (FIG. 7A) and of consecutive edges onface-adjacent hexes (FIG. 7B) are expected to be quite similar.Optimizing for direction similarity of such edges may help to guide thetarget edge directions across the mesh toward a more uniform solution,thereby improving the stability of the optimization process of methods100, 300. Some embodiments employ a regularization term of the form:

$\begin{matrix}{E_{regularize} = {{\sum\limits_{e \in E}{\frac{1}{{P(e)}}{\sum\limits_{e_{p} \in {P{(e)}}}{{\frac{e}{L_{e}} - \frac{e_{p}}{L_{p}}}}^{2}}}} + {\sum\limits_{e \in E}{\frac{1}{{C(e)}}{\sum\limits_{e_{c} \in {C{(e)}}}{{\frac{e}{L_{e}} - \frac{e_{p}}{L_{p}}}}^{2}}}}}} & (19)\end{matrix}$

where e_(p) ∈ P(e) are the edges topologically parallel to a given edgee and e_(c) ∈ C (e) are its consecutive edges. In the equation (19)formulation, all vectors are normalized by target edge lengths |L_(e)|,|L_(p)|, |L_(c)| to optimize for the desired length proportion.

A desirable feature in some mesh optimization applications and/orembodiments involves the preservation of the outer surface of the mesh(which may correspond to the outer surface of an object, for example).Some embodiments of methods 100, 300 may achieve this desirable featureby fixing the positions of the surface vertices. However, in someapplications and/or embodiments, such strict positional constraints aretoo restrictive, as often mesh quality can be improved by allowing thesurface vertices to slide along the boundary surface but without movingaway from the boundary surface. Moreover, in some applications and/orembodiments, one cannot obtain an inversion-free mesh without allowingthe mesh vertices to deviate from the boundary surface. Thus, in someembodiments, methods 100, 300 use a boundary formulation that balancessurface preservation against mesh quality.

To optimize the quality of hexahedra along the boundary surface, someembodiments include partial cones around each surface edge. As the aimmay be for surface edges to remain on the surface, the local step (e.g.block 360 of method 300) equation (9) may be modified for surface edges:

$\begin{matrix}{{\hat{n}}_{ij} = {{{argmin}{\overset{{T{(e_{ij})}}}{\sum\limits_{k = 1}}\left( {{\hat{n}}_{ij} - n_{k}} \right)^{2}}} + {\beta \; {n_{ij} \cdot \overset{\Cup}{n}}}}} & (20)\end{matrix}$

where: {hacek over (n)} is the average of the original surface normalsat v_(i) and v_(j); and β is a configurable constant (e.g. β=10). Thetreatment of boundary edges in the global step (e.g. step 370 of method300) may be similar to that of interior edges.

Expressing exact or approximate surface preservation in closed form(instead of just preserving the current locations) may involve the useof an analytic surface definition, which is typically rarely available.Thus to constrain vertices to remain on or close to the surface, someembodiments may employ a local surface approximation. Some embodimentsmay distinguish between regular surface vertices v ∈ S, feature verticesv ∈ F and corners (vertices at the junction of multiple features) v ∈ C.To preserve the surface, some embodiments may use three types ofper-vertex energy terms, selected by vertex type. Some embodiments mayuse a boundary energy term having the form:

$\begin{matrix}{E_{boundary} = {{\sum\limits_{v \in S}{\beta \left( {{\hat{n} \cdot v} + \hat{d}} \right)}^{2}} + {\sum\limits_{v \in F}\left( {{\alpha \left( {v - \left( {\hat{v} + {a\hat{t}}} \right)} \right)}^{2} + a^{2}} \right)} + {\sum\limits_{v \in C}{\alpha \left( {v - \hat{v}} \right)}^{2}}}} & (21)\end{matrix}$

where: {circumflex over (v)} is the reference, or closest, surfaceposition for each vertex; <{circumflex over (n)}, {circumflex over (d)}>is the implicit equation of the plane passing through {circumflex over(v)} and orthogonal to the input surface normal {circumflex over (n)} atthat point; {circumflex over (t)} is the feature tangent at {circumflexover (v)}; a is an auxiliary variable added to the system to enablefeature constraints; and α is a configurable constant for feature andcorner weight. In some embodiments, α is set to α=20. The equation (21)boundary term E_(boundary) may be added, for example, to the right handside of equation (1a), (1b) and/or (18). The equation (21) boundary termE_(boundary) may form part of the term Ê in equation (18). It may beobserved that the equation (21) boundary term E_(boundary) is quadratic,thus augmenting the combined energy function E_(global) with theequation (21) boundary term E_(boundary) does not affect the convexityof the global optimization performed in methods 100, 300.

As discussed above, methods 100, 300 may make use of target, or optimal,edge lengths in a number of places in the formulation described above.Such target edge lengths may be determined in block 120 of method 100and/or block 320 of method 300, for example. Such target edge lengthsmay be used in the remainder of method 100 and/or method 300. In manymeshing frameworks, such target edge lengths are provided by users andmay be based on simulation needs. Absent this high-level knowledge, whenthe mesh density is by design (or otherwise) non-uniform (e.g. formeshes created with octree-based methods), the input edge lengthsprovide a plausible estimation of the target edge length. If the meshdensity is expected to be uniform, some embodiments may provide a betterinitial length estimate by leveraging the existing surface mesh sizing.In general, one expects mesh size changes to be gradual, indicating apreference for similar lengths on both topologically parallel andconsecutive edges. Further, as a uniformly sized mesh is sought, theaverage edge length L_(A) in the current mesh (which may be the inputmesh) may serve as a good initial approximation of target interior meshedge lengths. Some embodiments involve combining these observations tomotivate an optimization having an edge-length energy functionE_(length) having the following formulation for obtaining target edgelengths L_(e):

E _(length=Σe∈Sr(L) _(e) _(−∥e∥) ²)²+Σ_(e∈E/Sr(L) _(e) _(−L) _(A) ₎ ₂_(+Σep∈p(e)(L) _(e) _(−L) _(p) ₎ ₂ +Σ _(e) _(c) _(∈c(e)(L) _(e) _(−L)_(c) ₎ ₂   (22)

where: e ∈ S′ are the surface edges having both vertices on a boundarysurface; e ∈ E/S′ are non-surface edges having at least one vertex thatis not on a boundary surface; e_(p) ∈ P(e) are the edges topologicallyparallel to e; e_(c) EC (e) are the edges topologically consecutive toe; ∥e∥ are original edge lengths; L_(A) is the average edge length inthe mesh; L_(p) is the length of a a topologically parallel edge e_(p);and L_(c) is the length of a topologically consecutive edge e_(c).Blocks 120 and/or 320 described above, may comprise optimizing (e.g.minimizing) an equation having the form of equation (22). The resultingedge lengths L_(e) balance the expectation of uniformity against thedesire to preserve the current surface edge lengths. The resulting setof edge lengths L_(e) determined in accordance with this edge-lengthoptimization may be used as the target edge lengths L_(ij) described atother locations herein.

In some embodiments, the overall hex-mesh optimization (e.g. of method100 and/or method 300) can start with these target edge-lengths (e.g. byestimating same in blocks 120 and/or 320). In some embodiments, thesetarget lengths can be updated selectively, reflecting the local meshelement quality, by linearly interpolating between the target edgelengths L_(e) and current edge lengths ∥e∥,

L _(e)=(1−w)·L _(e) +w∥e∥  (23)

where w is a Gaussian function of the smallest hex element minimalscaled jacobian (MSJ) q of any hex connected to the target edge:

$\begin{matrix}{w = ^{{- 0.5} \cdot {(\frac{q}{\sigma})}^{2}}} & (24)\end{matrix}$

where σ is a configurable constant (set, in some embodiments, toσ=0:15). This selective update aims to preserve the original targetlengths in areas where the mesh quality is sufficiently high, so as topreserve the mesh in these areas as is, while updating edge lengths inareas where the quality is poor.

As discussed above, in some embodiments, the global energy functionE_(global) optimized in method 100 or in the global procedure 370 ofmethod 300 may comprise a sum of some of the optional energy termsdiscussed above. In some embodiments, the global energy functionE_(global) optimized is given by equation (18), where the Ê term isgiven by:

Ê=E _(cone) +E _(regularize) +E _(boundary)   (25a)

where: E_(cone) has a form of equation (2) or equation (5), whereE_(regularize) has a form of equation (19) and E_(boundary) has a formof equation (21). This yields a global energy function E_(global) havingthe form:

E _(global) =E _(cone) +E _(penalty) +E _(regularize) +E _(boundary)  (25b)

where E_(penalty) is the summation in the right hand side of equation(18). Since the minimized functional (the global energy functionE_(global)) is quadratic, some embodiments make use of a standard linearsolver (GMRES or SuperLU, for example) to compute the minimizer ofequation (25). Once the iterations of the global procedure converge(e.g. iterative changes to the vertex positions v_(i) or the energyvalue E_(global) are less that a suitable threshold), some embodimentsmay optionally comprise repeating the hex-mesh optimization process(e.g. method 100 or method 300) with updated target edge lengths (asdiscussed above e.g. repeating method 100 and/or 300 with updates to thetarget edge lengths determined in optional block 120 and/or 320).

The framework described above is designed to obtain inversion-freemeshes with high average element quality. However, some simulations haveminimal, or worst, element quality restrictions that are more stringentthan mere convexity. In some embodiments, the formulation describedabove aims for minimal angle between each directed edge and the basetriangle normals of its associated edge-cone to be some configurableconstant θ close to 90° (e.g. θ=85°, in some embodiments) and usesε=cos(θ) in equation (10). In some embodiments, to achieve an improvedminimum quality, a variation of the above-described techniques may beused to progressively improve the minimal quality of the resultanthex-mesh. In particular, some embodiments may comprise changes toequations (10) and (16) to achieve this result. Once all local solutionspaces are not empty (i.e. all cone angles α_(ij) are below 90°, thedefault minimal angle threshold ε (in equation (10)) may be replacedwith a local threshold whose goal is to further improve the localcone-quality:

{circumflex over (n)} _(ij)·>cos(α_(ij)+ε) for all t _(k) ∈T(e _(ij))  (26)

where α_(ij) is the worst current angle between the directed edge e_(ij)and the normals of the base triangles of its corresponding edge-cone andε is some configurable constant (e.g. ε=0.01, in some embodiments). Thischange imposes a requirement that all directed edge target directions{circumflex over (n)}_(ij) determined by the local procedure (block 360)be better aligned with respect to the normals of the base triangles ofits corresponding edge-cone than they are using the above-describedequation (10). In some embodiments, the above-described target direction{circumflex over (n)}_(ij) may be kept if and when the equation (26)requirement results in an empty solution space. In some embodiments, noeffort is made to improve the local cones if cos(a_(u)) is above asuitably high-quality threshold (e.g. cos(α_(ij))>0.5).

Some embodiments also comprise determining the current worst angle wacross all edge-cones. If this current worst angle value w is under 90°,the angle-based weight W(α_(eij)) in the penalty function (equation(16)) may be updated to be maximized at this worst angle,

$\begin{matrix}{{W\left( \alpha_{eij} \right)} = {1 + {10^{\frac{- {({{\cos {(\alpha_{ij})}} - {\cos {(\omega)}}})}^{2}}{2\sigma^{2}}}}}} & (27)\end{matrix}$

Mesh-optimization loop 330 (with these updates for maximizing minimalquality (e.g. equations (26) and (27))) may be performed, re-determiningα_(ij) and w each time the process converges (e.g. for each iteration ofmesh-optimization loop 330) and may be concluded when furtherimprovement is no longer possible (e.g. w no longer increases) or oncethe minimal quality reaches a user-configurable minimal quality value.

FIG. 8 is a schematic depiction of a system 500 which may be used toperform any of the methods described herein and the steps of any of themethods described herein according to a particular embodiment. System500 of the illustrated embodiment comprises one or more computers 502which may comprise one or more processors 504 which may in turn executesuitable software (not expressly enumerated) accessible to processor(s)504. When such software is executed by computer 502 (and in particularprocessor(s) 504), computer 502 and/or processor(s) 504 may perform anyof the methods described herein and the steps of any of the methodsdescribed herein. In the illustrated embodiment, computer 502 provides auser interface 510 for interaction with a user 506. From a hardwareperspective, user interface 510 comprises one or more input devices 508by which user 506 can input information to computer 502 and one or moreoutput devices 512 by which information can be output to user 506. Ingeneral, input devices 508 and output devices 512 are not limited tothose shown in the illustrated embodiment of FIG. 6. In general, inputdevice 508 and output device 512 may comprise any suitable input and/oroutput devices suitable for interacting with computer 502. Userinterface 510 may also be provided in part by software when suchsoftware is executed by computer 502 and/or its processor(s) 504. In theillustrated embodiment, computer 502 is also connected to access data(and/or to store data) on accessible memory device 518. In theillustrated embodiment, computer 502 is also connected by communicationinterface 514 to a LAN and/or WAN network 516, to enable accessing datafrom networked devices (not shown) and/or communication of data tonetworked devices.

Input (e.g. input hex-mesh 102) may be obtained by computer 502 via anyof its input mechanisms, including, without limitation, by any inputdevice 508, from accessible memory 518, from network 516 or by any othersuitable input mechanism. The outputs (e.g. output hex-mesh 158) may beoutput from computer 502 via any of its output mechanisms, including,without limitation, by any output device 512, to accessible memory 518,to network 516 or to any other suitable output mechanism. As discussedabove, FIG. 8 is merely a schematic depiction of a particular embodimentof a computer-based system 500 suitable for implementing the methodsdescribed herein. Suitable systems are not limited to the particulartype shown in the schematic depiction of FIG. 8 and suitable components(e.g. input and output devices) are not limited to those shown in theschematic depiction of FIG. 8.

The methods described herein may be implemented by computers comprisingone or more processors and/or by one or more suitable processors, whichmay, in some embodiments, comprise components of suitable computersystems. By way of non-limiting example, such processors could comprisepart of a computer-based automated contract valuation system. Ingeneral, such processors may comprise any suitable processor, such as,for example, a suitably configured computer, microprocessor,microcontroller, digital signal processor, field-programmable gate array(FPGA), other type of programmable logic device, pluralities of theforegoing, combinations of the foregoing, and/or the like. Such aprocessor may have access to software which may be stored incomputer-readable memory accessible to the processor and/or incomputer-readable memory that is integral to the processor. Theprocessor may be configured to read and execute such softwareinstructions and, when executed by the processor, such software maycause the processor to implement some of the functionalities describedherein.

Certain implementations of the invention comprise computer processorswhich execute software instructions which cause the processors toperform a method of the invention. For example, one or more processorsin a computer system may implement data processing steps in the methodsdescribed herein by executing software instructions retrieved from aprogram memory accessible to the processors. The invention may also beprovided in the form of a program product. The program product maycomprise any medium which carries a set of computer-readable signalscomprising instructions which, when executed by a data processor, causethe data processor to execute a method of the invention. Programproducts according to the invention may be in any of a wide variety offorms. The program product may comprise, for example, physical(non-transitory) media such as magnetic data storage media includingfloppy diskettes, hard disk drives, optical data storage media includingCD ROMs, DVDs, electronic data storage media including ROMs, flash RAM,or the like. The instructions may be present on the program product inencrypted and/or compressed formats.

Where a component (e.g. a software module, controller, processor,assembly, device, component, circuit, etc.) is referred to above, unlessotherwise indicated, reference to that component (including a referenceto a “means”) should be interpreted as including as equivalents of thatcomponent any component which performs the function of the describedcomponent (i.e., that is functionally equivalent), including componentswhich are not structurally equivalent to the disclosed structure whichperforms the function in the illustrated exemplary embodiments of theinvention.

As will be apparent to those skilled in the art in the light of theforegoing disclosure, many alterations and modifications are possible inthe practice of this invention without departing from the spirit orscope thereof. For example:

-   -   This description uses the terms “global”, “globally” as a matter        of convenience to distinguish from operations that are “local”        or performed “locally”. Unless the context dictates otherwise,        the terms “global” or “globally” should not be interpreted        strictly to require that an operation performed globally is        performed for every vertex and/or for every edge of a mesh. For        example, in some cases, operations performed globally could        include operations performed on an entirety of current mesh 112        or input mesh 102 (e.g. for every vertex and/or edge of the        mesh), but this is not strictly necessary. In some embodiments,        or applications, operations performed “globally” may be        performed on a portion or portion(s) of current mesh 112 or        input mesh 102 that includes a plurality of directed edges or a        plurality of edge-cones. Similarly, in some embodiments, an        object represented by an input mesh 102, an input mesh 102        and/or a current mesh 112 could be parsed into several portions,        each of which could be treated as an input mesh 102 or a current        mesh 112 on which operations could be performed globally. By way        of non-limiting example, the methods described herein may be        applied to problem portions of the mesh, while leaving        non-problematic portions of the mesh unchanged.    -   The above-described penalty terms represent only one form of        penalty term that could be used in particular embodiments (e.g.        with energy function of equation (18) and/or equation (25)). In        other embodiments, other penalty terms/functions could be used.        Such penalty functions may express a preference for directed        edges being similar to the target directed edges determined in        the local procedure.    -   In some embodiments, block 120 or block 320 target edge lengths        may be selected to be current edge lengths. Such embodiments may        involve the use of a non-linear solver.

While a number of exemplary aspects and embodiments have been discussedabove, those of skill in the art will recognize certain modifications,permutations, additions and sub-combinations thereof. It is thereforeintended that the following appended aspects or claims and aspects orclaims hereafter introduced are interpreted to include all suchmodifications, permutations, additions and sub-combinations.

What is claimed is:
 1. A method, performed by a computer, for improvingquality of a hex-mesh, the method comprising: obtaining an inputhex-mesh; performing a first iterative procedure which starts with theinput hex-mesh and which outputs an output hex-mesh having improvedquality relative to the input hex-mesh; wherein performing the firstiterative procedure comprises: initializing a current hex-mesh to beequal to the input hex mesh; for each iteration of the first iterativeprocedure: performing an optimization of an energy function over aplurality of directed-edges in the current hex-mesh to determine updatedvertex positions for vertices in the current hex-mesh, wherein for eachdirected edge, the energy function comprises a term that expresses apreference for the directed edge to be aligned with normal vectors ofbase triangles of an edge-cone corresponding to the directed edge; andupdating the current hex-mesh with the updated vertex positions; andafter one or more iterations, setting the output hex-mesh to be equal tothe current hex-mesh.
 2. A method according to claim 1 wherein, for eachdirected edge extending from a base vertex v_(i) to an end vertex v_(j):the edge-cone corresponding to the directed edge comprises a pluralityof tetrahedra surrounding the directed edge, each of the plurality oftetrahedra bounded by the directed edge and by a pair of distinct edgesof a hexahedron in the current hex-mesh mesh, each of the pair ofdistinct edges emanating from the base vertex v_(i) and not includingthe end vertex v_(j); the base triangles of the edge-cone correspondingto the directed edge comprise, for each tetrahedron of the plurality oftetrahedra, the triangular surface of the tetrahedron which includes thebase vertex v_(i) but does not include the end vertex v_(i).
 3. A methodaccording to claim 1 wherein, for each directed edge, the term of theenergy function that expresses the preference for the directed edge tobe aligned with normal vectors of base triangles of the edge-conecorresponding to the directed edge has the form of${\sum\limits_{t_{k} \in {T{(e_{ij})}}}\left( {\frac{e_{ij}}{e_{ij}} - n_{k}} \right)^{2}},$as described above in connection with equation (2).
 4. A methodaccording to claim 1 comprising, for each directed edge, at least oneof: determining the normal vectors of the base triangles correspondingto the directed edge according to an equation having the form ofequation (4) described above; and obtaining target edge lengths for theoutput hex-mesh and determining the normal vectors of the base trianglescorresponding to the directed edge according to an equation having theform of equation (6) described above.
 5. A method according to claim 4wherein obtaining target edge lengths for the output hex-mesh comprisesone or more of: receiving the target edge lengths for the outputhex-mesh; using the edge-lengths of the input hex-mesh as the targetedge lengths for the output hex-mesh; and estimating the target edgelengths for the output hex-mesh.
 6. A method according to claim 5wherein obtaining target edge lengths for the output hex-mesh comprisesestimating the target edge lengths for the output hex-mesh andestimating the target edge lengths for the output hex-mesh comprisesperforming an edge-length optimization which optimizes an edge-lengthenergy function over a plurality of edges in either the input hex-meshor the current hex-mesh.
 7. A method according to claim 6 where theedge-length energy function has the form of equation (22) describedherein.
 8. A method according to claim 1 wherein performing theoptimization of the energy function comprises performing an iterativeoptimization and wherein, for each iteration of the iterativeoptimization, the energy function comprises, for each directed edge, aweighted penalty term which expresses a preference for the directed edgeto be aligned with a target edge direction {circumflex over (n)}_(ij).9. A method according to claim 8 wherein, for each iteration of theiterative optimization, updating the weighted penalty term for eachdirected edge.
 10. A method according to claim 8 wherein performing thefirst iterative procedure comprises, for each iteration of the firstiterative procedure and for each directed edge, determining the targetedge direction {circumflex over (n)}_(ij).
 11. A method according toclaim 10 wherein, for each directed edge, determining the target edgedirection {circumflex over (n)}_(ij) is based on characteristics of thecurrent hex-mesh local to the edge-cone corresponding to the directededge.
 12. A method according to claim 11 wherein, for each directededge, the characteristics of the current hex-mesh local to the edge-conecorresponding to the directed edge are the normal vectors of the basetriangles of the edge-cone corresponding to the directed edge and thepositions of the vertices which define the directed edge.
 13. A methodaccording to claim 10 wherein, for each directed edge, determining thetarget edge direction {circumflex over (n)}_(ij) comprises performing alocal optimization which optimizes a local energy function comprising aplurality of local terms, each local term expressing a preference forthe target edge direction {circumflex over (n)}_(ij) to be aligned witha normal vector of a corresponding one of the base triangles of theedge-cone corresponding to the directed edge.
 14. A method according toclaim 10 wherein, for each directed edge, determining the target edgedirection {circumflex over (n)}_(ij) comprises: initially setting aninitial target edge direction {circumflex over (n)}_(ij) to correspondto the direction of an average of the normal vectors of the basetriangles of the edge-cone corresponding to the directed edge;evaluating whether the initial target edge direction {circumflex over(n)}_(ij) satisfies non-inversion constraints; and if the initial targetedge direction {circumflex over (n)}_(ij) satisfies non-inversionconstraints, then determining the target edge direction {circumflex over(n)}_(ij) to be the initial target edge direction {circumflex over(n)}_(ij) and, otherwise, performing a local optimization whichminimizes a local energy function which comprising a plurality of localterms, each local term expressing a preference for the target edgedirection {circumflex over (n)}_(ij) to be aligned with a normal vectorof a corresponding one of the base triangles of the edge-conecorresponding to the directed edge.
 15. A method according to claim 10wherein, for each iteration of the iterative optimization, the weightedpenalty term is based on the target edge direction {circumflex over(n)}_(ij).
 16. A method according to claim 15 wherein, for each directededge, the weighted penalty term comprises a penalty function E_(ij) ^(p)which is optimized as part of the energy function and a penalty weightw_(ij) ^(p) that is updated between iterations of the second iterativeoptimization.
 17. A method according to claim 16 wherein, for eachdirected edge, the penalty function E_(ij) ^(p) expresses the preferencefor the directed edge to be aligned with a target edge direction ñ_(ij).18. A method according to claim 15 wherein, for each directed edge, theweighted penalty term comprises a geometry factor W(α_(ij)) based atleast in part on a largest angle α_(ij) between the target direction{circumflex over (n)}_(ij) and the normal vectors of base triangles ofthe edge-cone corresponding to the directed edge.
 19. A method accordingto claim 8 wherein, for each iteration of the iterative optimization,the energy function comprises a parallel-edges regularization termE_(regularize) _(_) _(parallel), which expresses a preference fordirections of topologically parallel edges of the same hexahedron to bethe same.
 20. A method according to claim 8 wherein, for each iterationof the iterative optimization, the energy function comprises aconsecutive-edges regularization term E_(regularize) _(_)_(consecutive), which expresses a preference for directions oftopologically consecutive edges of face-adjacent hexahedra to be thesame.
 21. A method according to claim 13 wherein, for each surfacedirected edge having both of its vertices on a surface boundary of theinput mesh, performing the local optimization comprises optimizing aboundary-edge energy function which expresses a preference for thesurface directed edge to remain on the surface boundary.
 22. A methodaccording to claim 8 wherein the energy function comprises a boundarypreservation term E_(boundary) which expresses a preference for verticeswhich are on a surface boundary of the input hex mesh to remain on theboundary surface.
 23. A method according to claim 13 wherein, for eachdirected edge, performing the location optimization comprises performingthe local optimization with first non-inversion constraints until thecurrent mesh is free from inverted hexahedra and then switching thefirst non-inversion constraints to first quality-improvement constraintswhich require improved quality relative to previous iterations of thelocal optimization.
 24. A method according to claim 18 comprising, foreach iteration of the first iterative process, determining a worst anglew from among the largest angles α_(ij) across the plurality ofedge-cones corresponding to the plurality of directed edges in thecurrent mesh, and if w is less than 90°, changing the geometry factorW(α_(ij)) from a first geometry factor to a second geometry factor forsubsequent iterations of the first iterative process, the secondgeometry factor different from the first geometry factor.
 25. A systemfor improving quality of a hex-mesh, the system comprising: one or moreprocessors configured to: obtain an input hex-mesh; perform a firstiterative procedure which starts with the input hex-mesh and whichoutputs an output hex-mesh having improved quality relative to the inputhex-mesh; wherein the one or more processors are configured to performthe first iterative procedure by: initializing a current hex-mesh to beequal to the input hex mesh; for each iteration of the first iterativeprocedure: performing an optimization (e g minimization) of an energyfunction over a plurality of directed-edges in the current hex-mesh todetermine updated vertex positions for vertices in the current hex-mesh,wherein for each directed edge, the energy function comprises a termthat expresses a preference for the directed edge to be aligned withnormal vectors of base triangles of an edge-cone corresponding to thedirected edge; and updating the current hex-mesh with the updated vertexpositions; and after one or more iterations, setting the output hex-meshto be equal to the current hex-mesh.